The mechanical response of semifiexible networks to localized perturbations 
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Previous research on semifiexible polymers including cytoskeletal networks in cells has suggested 
the existence of distinct regimes of elastic response, in which the strain field is either uniform 
(affine) or non-uniform (non-affine) under external stress. Associated with these regimes, it has 
been further suggested that a new fundamental length scale emerges, which characterizes the scale 
for the crossover from non-affine to affine deformations. Here, we extend these studies by probing 
the response to localized forces and force dipoles. We show that the previously identified nonaffinity 
length [D.A. Head et al. PRE 68, 061907 (2003).] controls the mesoscopic response to point forces 
and the crossover to continuum elastic behavior at large distances. 

PACS numbers: 87.16.Ka, 82.35.Pq, 62.20.Dc 



I. INTRODUCTION 

Semifiexible polymers such as filamentous proteins re- 
semble elastic rods on a molecular scale, while exhibiting 
significant thermal fluctuations on the scale of microm- 
eters or even less. This has made them useful as model 
systems allowing for direct visualization via optical mi- 
croscopy. But, semifiexible polymers are not just large 
versions of their more well-studied flexible cousins such 
as polystyrene. Filamentous proteins, in particular, have 
been shown to exhibit qualitatively different behavior in 
their networks and solutions. A fundamental reason for 
this is the fact that the thermal persistence length, a mea- 
sure of filament stiffness as the length at which thermal 
bending fluctuations become apparent, can become large 
compared with other important length scales such as the 
spacing between polymers in solutions, or the distance 
between chemical crosslinks in a network. 

One of the most studied semifiexible polymers in re- 
cent years has been F-actin, a filamentous protein that 
plays a key structural role in cells Q, 0, Q • These oc- 
cur in combination with a wide range of specific proteins 
for crosslinking, bundling and force generation in cells. 
These composites, together with other filamentous pro- 
teins such as microtubules, constitute the so-called cy- 
toskeleton that gives cells both mechanical integrity and 
structure. This biopolymer gel is but one example of a 
large class of polymeric materials that can store elastic 
energy in a combination of bending and extensional de- 
formations of the constituent elements. Such systems can 
be called semifiexible gels or networks. 

One of the important lessons from recent experimen- 
tal and theoretical studies is that the shear modulus 
of cross-linked semifiexible networks bears a much more 
complex relationship to the mechanical properties of the 
constituent filaments and to the microstructure of the gel 
than is the case for flexible polymer gels . Recently, it 
has been shown that semifiexible gels exhibit a striking 



cross-over [E IS H, between a response to external 
shear stress that is characterized by a , spa tially heteroge- 
neous strain (a non-affine regime llfjlllll ) and a uniform 
strain response (an affine regime |l2|). This crossover 
is governed primarily by cross-link density and molecu- 
lar weight (filament length). The bulk shear modulus 
of the network simultaneously increases by about six or- 
ders of magnitude at this nonaffine-to-affine cross-over. 
The underlying mechanism responsible for this abrupt 
cross-over appears to be the introduction of a new, meso- 
scopic length scale in the problem that is related to both 
the bending stiffness of the constituent polymers and the 
mean spacing between consecutive cross-links along the 
chain [EH- 

One can associate this mesoscopic length with the 
length below which the deformation of the network de- 
parts from the standard affinity. The nonaffinity length 
A, introduced in Refs. can be qualitatively under- 

stood as the typical length over which one finds nonaffme 
deformation in the network. In this previous work we 
presented a scaling analysis that relates this mesoscopic 
length scale to the network density and the stretching and 
bending moduli of the constituent filaments. The macro- 
scopic shear response of the network is then controlled 
by a competition between A (the nonaffinity length) and 
the filament length, L. On the one hand, when the fil- 
ament length is long, nonaffme corrections to the de- 
formation field, which are localized to regions within A 
of the filament ends, do not significantly affect the me- 
chanical properties of the network; the shear modulus of 
the macroscopic system is well-described by calculations 
based on affine deformation reflecting the fact that non- 
affine deformations of the filament ends are subdominant 
corrections in this limit. Moreover, the elastic energy 
is stored primarily in the (homogeneous) extension and 
compression of filaments. On the other hand, when the 
filaments are of a length comparable to, or shorter than 
the nonaffinity length, i.e. L ~ A then the nonaffme de- 
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formations of the ends play a large, even dominant role 
in determining the mechanical properties of the network. 
The network is found to be generally more compliant, and 
the elastic energy under applied shear stress is stored in 
a spatially heterogeneous manner in the bending of fila- 
ments. The existence of these distinct regimes as a func- 
tion of filament length reflects a fundamental difference 
of these semiflexible polymer networks with respect to 
their flexible counterparts: polymers can maintain their 
mechanical integrity and state of stress/strain across net- 
work nodes or crosslinks. 

These results naturally lead one to pose a number of 
basic questions regarding the elastic properties of semi- 
flexible networks. While these random networks must on 
basic theoretical grounds appear as continuum, isotropic 
materials at the longest length scales, these considera- 
tions do not predict the length at which the continuum 
approximation applies. The previous identification of the 
nonaffinity length, A as the only mesoscopic length associ- 
ated with the nonaffine-to-affine cross-over in uniformly 
sheared semiflexible networks suggests that this length 
should more generally control the cross-over to contin- 
uum behavior ||| . After all the affine deformation of the 
network under uniform stress at scales large compared to 
A shows that in one case at least the nonaffinity length 
controls the cross-over to continuum behavior. One of 
the principal results of the present work is the demon- 
stration that A more generally controls this cross-over to 
continuum mechanics in semiflexible gel systems. 

Prior work has focussed exclusively on simple shear 
and uniaxial extension. In order to better examine the 
universality of the previous results, we study the oppo- 
site limit of a highly localized external force in the form 
of a point force monopole or dipolc. If one can show 
that the elastic (displacement) Green's function of the 
system similarly depends on only one additional param- 
eter A then it would appear that this quantity completes 
the coarse-grained elastic description of the system on all 
length scales down to the mean distance between cross- 
links. It may be, however, that the deformation field of 
these semiflexible networks is much more complex and 
the simplification introduced by A in the description of 
the network's response to uniform shear strain cannot be 
generalized to deformations resulting from more general 
forcing conditions. 

While our results below, indeed, show that A docs 
largely control the crossover to (the far field) contin- 
uum elasticity, the observed elastic Green's function is 
sensitive to the local structure of the network on length 
scales below A. Below we discuss how we quantify the 
structure of the Green's function and its approach to the 
form required by continuum elasticity. The observed elas- 
tic Green's function, however, depends not only on the 
A-dependent Lame coefficients of the material, but also 
on local properties of the displacement field immediately 
surrounding the the point force. In effect one may imag- 
ine that, upon the application of the point force, the 
network acts as a type of composite material: within a 



distance A of the point force it deforms in a way not well 
described by continuum theories, while outside of that 
zone it does appear to act like an elastic continuum. The 
complete Green's function depends, of course, on the ma- 
terial properties of both media. Unfortunately, only one 
of those media (the outer zone) is well characterized by 
the simple continuum theory of an isotropic elastic solid, 
so the complete Green's function remains complex and 
depends on the detailed network structure within the in- 
ner zone surrounding the applied force. 

There is another set of questions that may be addressed 
via the study of the network's response to point forces 
and force dipoles. Such forces not only probe the mate- 
rial properties of the network in a manner complimentary 
to the uniform strains explored earlier, they also have di- 
rect physical implications for microrheology in F-actin 
networks and for the dynamics of the cytoskeleton in re- 
sponse to the activity of nanoscale molecular motors, e.g. 
myosin. Fluctuation-based microrheology, an application 
of the fluctuation-dissipation theorem to the study of rhe- 
ology via the statistical analysis of the thermally fluctu- 
ating position of sub-micron tracer particles embedded 
in the medium, requires one to understand the elastic 
Green's function of the medium. Thus understanding 
the response of semiflexible networks to localized forces 
has direct experimental implications and consequences 
for force production in the cytoskeleton. 

In the biological context, the semiflexible network mak- 
ing up the cytoskeleton is generally found in association 
with molecular motors that, to a good approximation, 
generate transient localized force dipoles in the material. 
To both understand force generation in the cell as well 
as the material properties of these cytoskeletal networks 
driven into nonequilibrium steady states by these molec- 
ular motors, one must determine the displacement field 
associated with such motor-induced forces. 

Notwithstanding our biological motivation for this 
work, our findings also bear on the broader problem of 
elastic modes in amorphous materials. It has been shown 
that the vibrational modes of deep-quenched Lennard- 
Jones systems approach a continuum description only on 
scales exceeding some mesoscopic length £ 0, Q ; for 
the protocols considered, a value £ ~ 30 particle dimen- 
sions was robustly found. This was physically identified 
with a length scale for non-affinity, suggesting a direct 
correspondence with our A (although our A can be con- 
trolled by varying the mechanical properties of the con- 
stituents) . A comparable length was also found to control 
the self-averaging of the Green's function to the form ex- 
pected by continuum elasticity (lE) . These findings for 
radially interacting particles are broadly in keeping with 
our own investigations for semiflexible polymer networks. 
We also mention here that the relationship between con- 
tinuum elasticity and the Green's function has also been 
discussed for mildly disordered spring networks 0] . 

The remainder of this paper is organized as follows: In 
section [n] we develop our model of semiflexible, perma- 
nently cross-linked gels, summarize the numerical Simula- 
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tions used to study it, and discuss the expected structure 
of the displacement field when averaged over numerous 
realizations of the network. In section ITTT1 we report our 
results for both point forces in IIII Al and force dipoles 
in IIII Bl We then discuss our studies of the bulk elastic 
properties of these networks in llll Cl bcfore concluding in 
section IIVI 

II. MODEL 
A. The semiflexible network 

A highly successful continuum continuum model of in- 
dividual semiflexible polymers is the worm-like chain. 
This treats the linear filaments as elastic rods of fixed 
contour length and negligible thickness, so that the dom- 
inant contribution to the elastic energy comes from bend- 
ing modes and the Hamiltonian linearized to small devi- 
ations from a straight configuration that is given by 

= \n J{V 2 ufds, (1) 

where u is the transverse displacement of the filament 
relative to an arbitrary straight axis, s is its arc length, 
and the elastic modulus k gives the bending energy per 
unit length Ss. 

The longitudinal response of the wormlike chain model 
is calculated from the increase in free energy due to an 
extensional stress applied along the mean filament axis 
|l2l Il7|. However, the numerical algorithm used in our 
simulations is based on the minimization of the Hamilto- 
nian of the system, and hence is fundamentally athermal. 
The reason for this choice is essentially one of efficiency: 
assuming there are no subtle convergence issues, mini- 
mization is expected to be faster than stochastic mod- 
elling and, it is anticipated, give better statistics for a 
given CPU time. This does, however, mean that the en- 
tropic mechanism governing the longitudinal response is 
absent, and an explicit energetic term is required. 

An unconstrained filament at T = forms a straight 
configuration, and thus elongation of its end-to-end dis- 
tance must be accompanied by a change in the absolute 
contour length. It is therefore natural to incorporate lon- 
gitudinal modes by adding a second elastic term to the 
Hamiltonian for the extension or shortening of the fila- 
ment backbone, 

where dl(s)/ds gives the strain or relative change in lo- 
cal contour length, and fi is the Young's modulus of 
the filament (essentially a spring constant normalized to 
1/ [length]). Of course, such modes also exist in thermal 
systems, but may be dominated by the entropic spring 
terms |l2j . except possibly for very short filament seg- 
ments or very densely cross-linked gels. The connec- 



tion between thermal/entropic and athermal longitudinal 
compliance is discussed in greater detail in Ref. [8j- 

The two elastic c oeffic ients k and /i together define a 
length scale = \J [i, which will shall refer to as the 
intrinsic bending length by observing that an isolated fil- 
ament constrained to have different tangents at its end 
points will deform with this characteristic length. To 
avoid potential confusion, however, we note that this is 
not the typical length scale for bending deformations of 
a semiflexible filament. Rather, the bending energy of 
filaments tends to make the longest unconstrained wave- 
length bending mode the dominant one. Thus, for in- 
stance, in a crosslinked gel, filaments are expected to be 
bent primarily on a length comparable to the distance 
between crosslinks. Nonetheless lh is a useful measure 
of filament rigidity, in that large lh corresponds to rigid 
filaments, and small lh to flexible ones. 

Although k and /i have been introduced as fundamen- 
tal coefficients, if the filament is regarded as a continuous 
elastic body with uniform cross section at zero temper- 
ature, then they can both be expressed in terms of the 
characteristic filament radius a and intrinsic bulk mod- 
ulus if as K ~ Yfa 4 and [i ~ Yfa 2 . Thus lh ~ a, and 
thinner filaments are more flexible than thick ones (as 
measured by It,), as intuitively expected. Given the pos- 
sibility of entropic effects in fj,, however, we shall treat 
these as independent parameters of the theory y| 

The gel is constructed by depositing filaments of 
monodisperse length L and zero thickness onto a two- 
dimensional substrate. The center-of-mass position vec- 
tor and orientation of the filaments are uniformly dis- 
tributed over the maximum allowed range, so the system 
is macroscopically isotropic and homogeneous. When- 
ever two filaments overlap they are cross-linked at that 
point. Deposition continues until the required mass den- 
sity, as measured by the mean distance between crosslinks 
l c , has been reached. The network thus constructed can 
be described by three lengths: L, lh, and l c and one mod- 
ulus scale: \i. In two-dimensions l c characterizes both 
the mass density and cross-link density in spatially ran- 
dom, isotropic networks. The length lh characterizes the 
mechanical properties of the constituent filaments via the 
ratio of their bending to stretching compliance. The over- 
all modulus scale /i will be absorbed into the point forces 
applied to the network. 

Previously [H IU, we identified an additional length 
A = l c (l c /lb) z j where z ~ 1/3. This non-affinity length 
characterized the crossover from non-affme to affine net- 
work response. Specifically, for filament lengths L much 
larger than A {i.e., high molecular weight), the bulk net- 
work properties could be understood quantitatively in 
terms of affine strains, while significant non-affine effects 
were observed for L < A. Although this length arose nat- 
urally from considerations of bulk network properties [j| , 
its dependence on l c and lh can be understood in terms 
of a balance of stretching/compression and bending ener- 
gies of a single filament, treated in a self consistent way 
within a network |8(. It is important to note that this 
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(material) length is intermediate, between the (geomet- 
ric) network length l c and the macroscopic scale. In fact, 
in dilute networks, for which l c ^> lb, we shall argue be- 
low that the network can be thought of on scales > A as 
a quasi-continuum: continuous, as opposed to discrete, 
but not necessarily described by macroscopic continuum 
elasticity. We shall find that this non-affinity length will 
play a key role in our analysis of the displacement field 
on this intermediate/quasi-continuum scale. 

In order to apply the point forces, a crosslink is chosen 
at random and identified as the origin of the system; it is 
this crosslink that will later be perturbed. A fixed circu- 
lar boundary at radius R from the origin is imposed, and 
any filaments or filament segments that extend beyond 
the boundary are simply removed or truncated, respec- 
tively. Filaments ending on the rigid boundary are fixed 
there by another freely rotating bond as are found at 
all cross-links in the system. The allowed free rotation 
at the boundary means that the boundary supplies arbi- 
trary constraint forces on the network but cannot support 
any localized torques. 



B. Numerical method 

Details of the simulation method have been presented 
elsewhere jSj- Here we briefly summarize the procedure, 
with particular attention on those aspects that are cen- 
tral to the problems studied in this paper. 

The system Hamiltonian 7i({xj}) is constructed from 
discrete versions of and applied to the geometry 
generated by the random deposition procedure described 
above. The degrees of freedom {xi}, or 'nodes,' are the 
position vectors of both crosslinks and midpoints between 
crosslinks, the latter so as to incorporate the first bend- 
ing mode along the filament segments. Additional nodes 
could be included at the cost of additional run time, but 
are expected to have a small effect, since the longest un- 
constrained wavelengths tend to dominate bending de- 
formations. Crosslinks are treated as constraints on the 
relative position of each connected filament segment, but 
not on their relative rotation. Physically this corresponds 
to an incxtensible but freely rotating linkage. As previ- 
ously noted, constrained bending at crosslinks has a small 
effect except at high network concentrations (specifically, 
when l c becomes comparable to a) Nodes on the 
boundary are immobile. Note that, as in our earlier 
work 0, @j the network is assumed to be initially 
unstressed on both macroscopic and microscopic length 
scales. 

There are two ways in which the system may be per- 
turbed. The first, which we call monopole forcing, is 
to apply an arbitrarily small external force c)f ext to the 
crosslink at the origin. The network is then allowed to 
relax to a new configuration consistent with mechanical 
force balance at every node. The second, which we call 
dipole forcing, is to introduce a geometrical defect into 
the system by moving the central crosslink along the con- 



tour length of one of the filaments to which it belongs, 
but not the other. Physically, this corresponds, e.g., to 
a motor introducing relative motion of one filament with 
respect to another filament. In the simulations, this ef- 
fect is incorporated by infinitesimally 'shifting' the image 
of the central node with respect to other nodes on one 
filament. 

Once the perturbation has been specified, the dis- 
placements of the nodes in the new mechanical equilib- 
rium are calculated by minimizing the system Hamilto- 
nian 7i({xi}) using the conjugate gradient method [l8| . 
This generates a displacement field for the particular ge- 
ometry under consideration, as shown in Fig. ^ Note 
that 7i({xi}) is linearized about small nodal displace- 
ments {#Xi} from their original positions {x^}, so linear 
response is assured. The bulk response to non-linear 
strains has been recently studied by Onck et al. jl9j| . 




FIG. 1: (Color online) An example of a network of filaments 
of uniform length L (grey line segments) perturbed by an ex- 
ternal force, denoted by the large red arrow, applied to the 
crosslink at the origin of a circular system of radius R — L. 
The arrow lengths are logarithmically calibrated to the mag- 
nitude of the displacement of at each crosslink. In this exam- 
ple, L/l c w 29.1, X/L ~ 0.191 and the force is perpendicular 
to one of the filaments that form the central crosslink; forces 
can also be applied parallel to a filament. 



C. Decomposition of mean displacement field 

As shown in Fig.Q] the displacement field for a particu- 
lar network is quite complex and generically shows anti- 
correlations between displacements and the local mass 
distribution. Although these fluctuations reflect inher- 
ent and possibly interesting physical properties of the 
gels, a more basic and immediately applicable quantity 
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to measure is the mean displacement field, found by aver- 
aging many individual runs with differing geometries but 
identical system parameters i?, L, l c and A (or equiva- 
lently R, L, l c and Two examples for differing A (the 
radius of the green circle centered on the point of force 
application) are given in Fig.0 These plots demonstrate 
one of the primary results of this paper, namely that the 
crossover between continuum (or quasi-continuum) re- 
sponse at large lengths, to a more exotic displacement 
field at shorter lengths, happens at a length of order (A) 
with a prefactor close to unity. 

In order to precisely describe the structure of the defor- 
mation field, we consider its most general possible form. 
For monopole forcing, the displacement field u(r) at po- 
sition vector r relative to the origin can be projected onto 
3 other vectors, namely the direction of the external force 
f , the unit position vector f and the axis n of one of the 
filaments to which the crosslink is attached, 

ui = Gijfj , 
Gn = -L {g {r) r t f 3 + g^hihj + gWsA (3) 

where / is the magnitude of the external force, and /i a g 
is the shear modulus as predicted for affine deformation. 
This depends on L and l c but not A, and is included here 
to factor out the density dependence of the response. As 
defined, the g('' (and the below) are dimensionless 
quantities in two dimensions. 

Each g('' can be further decomposed into angular 
modes in 9, where cos# = n • f , 

gU=gP + 2 ]T g^cos(md). (4) 

m>0, 

Terms in sin(m#) vanish since the ensemble-averaged re- 
sponse must be invariant under 9 <-* —9, and cos(m#) 
terms with m odd also vanish due to n «-> — n invari- 
ance. This latter symmetry may appear to violate the 
known polarity of typical semiflexible biopolymers such 
as F-actin and microtubules pflj : however, we are inter- 
ested here in the mechanical properties of the filaments, 
which, within the approximations of our model, are in- 
deed invariant under n <-> — n. 

For dipole forcing, the above procedure is followed with 
the additional simplification that n is proportional to f , 
since the displacement (and hence force) dipole induced 
by the motion of a motor will always be parallel to one 
filament axis. The decomposition is therefore somewhat 
simpler, 

Hij = -LlhWhtj + hfriSij} . (5) 

The scalar / is the magnitude of the force dipole; since it 
is actually a displacement that is imposed, / is unknown 
as will be treated as a fitting parameter. The angular 




(b) 

FIG. 2: (Color online) Mean response after averaging O(10 5 ) 
networks with the same parameters as in Fig. Q (including 
one filament fixed perpendicular to the external force) except 
for A, which is (a) X/L « 0.089 and (b) X/L w 0.42. For easy 
visualization, a green circle of radius A has been inserted into 
the background of each plot. Vectors near the center of each 
system have not been plotted for clarity. 



decomposition is identical to before, 

= h { ] + 2^4' cos(m<9) . (6) 

m>0, 

Later sections will refer to the continuum solution for 
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each of the two forms of forcing. These are given in the 
appendix. For the monopole case, only two continuum 

(r) ( f) 

modes are non-zero, namely 5^ and g$ ■ We shall refer 
to these components of the Green's function as contin- 
uum modes in order to distinguish from those compo- 
nents {non-continuum modes) that must vanish in the 
continuum. We do this even though for our finite sys- 
tems the non-continuum modes do not vanish. 

The dipole modes are slightly more subtle: after aver- 
aging of many dipole fields generated by the simulation, 
the resulting field is quadrupolar. This is an immediate 
consequence of the means of forcing the system. Recall 
that in an elementary step, a motor moves parallel to 
one filament axis. This has the effect of compressing the 
filament in front of the dipole, while stretching the trail- 
ing segment. Thus two filament segments are perturbed, 
each of which can be treated as a force dipole which, for 
a particular network geometry, will be of different mag- 
nitudes and hence the resulting field is dipole. However, 
the net bias of this dipole is symmetrically distributed 
around zero, and thus vanishes after averaging, leaving 
a quadrupole field as shown in Fig. [3] As derived in the 

(r) (r) 

appendix, the non-zero modes for this field are h \ h 2 , 
h ( f) and h 2 f) . 




FIG. 3: (Color online) Mean displacement field after averag- 
ing over O(10 J ) individual fields induced by imposing a dis- 
placement dipole at the origin. The orientation of the dipole 
is given by the red arrow. As explained in the text, the mean 
dipole moment vanishes after averaging and the resulting field 
is quadrupole. A green circle at a radius A = 0.191L from the 
origin is also shown. 



III. RESULTS 

The mechanical response to a localized perturbation 
depends on the distance from the point of perturbation. 
We divide the discussion of these results into the follow- 
ing parts. First, we examine the response to a point 
force using the decomposition of the displacement field 
outlined in previous section. We contrast the decay of 
the non-continuum modes of the displacement field with 
the behavior of the continuum modes of the displacement 
field and discuss the finite-size effects in the simulations. 
At large length scales, we find that the deformation field 
approaches a quasi-continuum form, in which all but a 
small number of ensemble-average modes decay rapidly 
toward zero. The remaining modes of the deformation 
field are the same as those predicted by continuum elas- 
ticity. In other words, the strain field about a point force 
reflects the expected tensorial character and rotational 
symmetries based on continuum elasticity theory. We 
find that A again plays a central role in controlling the 
cross-over from the near-field to the quasi-continuum. 

We then turn to the spatial structure of the elastic en- 
ergy density field around the point force paying partic- 
ular attention to the partitioning of that energy density 
between stretching and bending modes of the filaments. 
We observe that the ratio of these energy contributions 
achieves the bulk value over much shorter distances from 
the point force than does the structure of the displace- 
ment field acquire its far-field, or bulk structure. We 
then extend our analysis to consider force dipoles in the 
medium. Lastly, we extend our previous analysis of the 
bulk elasticity of the filament network by examining both 
the Young's modulus and the Poisson ratio of the net- 
work. 



A. The response to point forces 

1. The Displacement Field 

In this section we focus on the short length scale be- 
havior of the monopole response, for which the non- 
continuum modes are non-zero. We wish to distin- 
guish two distinct forms of this non-continuum behav- 
ior: (i) higher angular modes gm) and gm with m > 
are non-zero, and (it) the gm modes do not vanish, i.e. 
the response depends on the orientation of the filament 
to which the force is applied. This latter observation 
gives a clear indication of how the response can 'see' the 
microscopic structure of the gel on short length scales. 

An example demonstrating the appearance of non- 
continuum modes at short lengths is given in Fig. 0] 
( f) 

which shows the g 2 mode for systems with different 
crosslink densities L/l c but with the filament flexibility 
chosen to give the same A w 0.191L in each case. This 
plot shows the decay of the cos(26>) amplitude of the com- 
ponent of the displacement field in the direction of the 
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applied force (at the origin). In all three systems one ob- 
serves a rapid decay of this angular harmonic that must 
vanish for a continuum isotropic system. Moreover, the 
characteristic length scale for this decay appears to be of 
order A (~ 0.2R), although we examine this point more 
quantitatively below. For all network densities and values 
of A studied it is clear that the magnitude of g 2 van- 
ishes rapidly with distance from the point force. No data 
are shown for larger values of L/l c since at high network 
densities the numerical convergence of the strain field is 
so slow as to prevent attaining meaningful statistics. 
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FIG. 4: g% (which is dimensionless) versus distance from 
point of force application r for systems with the same X/L ~ 
0.191 and crosslink densities L/l c as given in the key (L is the 
filament length). The system radius R — L in all cases. 

This picture of non-continuum modes decaying to near 
zero at a length comparable to A holds also for the m > 2 
modes of gil and for all the m > modes of g-m ■ For 
reasons of space we do not show these data here. Instead 
we present data in Fig. that demonstrate the decay 

(n) 

of (?g with r for three networks densities such that in 
each case A w 0.191Z. This non-continuum mode of the 
displacement field measures the circularly averaged am- 
plitude of the component of displacement in the direction 
n, i.e. along the axis of the rod to which the force has 
been applied. Clearly, such behavior has no counterpart 
in an isotropic continuum clastic model. It is interesting 
to note that this amplitude also appears to decay expo- 
nentially with a decay length of order A. 

In order to test quantitatively whether A indeed con- 
trols the decay of those components of the displacement 
field that have no counterparts in the continuum theory, 

( f) 

we examine, as an example, g 2 vs. r more closely in 
Fig. [S] Here we plot this for a range of values of net- 
works density and of A. If, as suggested above, the decay 
is exponential with characteristic length A, then plotting 
these data log-linear with radial distances scaled by A 
should cause all these curves to exhibit the same slope. 
We have shifted the data sets vertically in order to fa- 
cilitate visual comparison. (We note that, although a 
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FIG. 5: g { n) versus r/R for X/L « 0.191, R = L and the L/l c 
given in the key. 



known force is applied to the origin of our sample, since 
we expect the unknown local constitutive relations in our 
system to depend on density and other parameters, the 
amplitudes cannot be directly compared. Here, we wish 
only to establish the nature of the decay of the non- 
continuum modes.) We have also introduced the solid 
lines corresponding to exp(— r/X) merely as guides to 
the eye. These are not fits, although fits to these data 
for various L/l c and A demonstrate decay lengths of A 
within 10%. 
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FIG. 6: Log linear plots of g 2 versus the scaled distance 
from perturbation r/X for systems with various values of A as 
given in the key. The system radius R = L in all cases. Data 
within a boundary layer of width A from r — R is not shown 
(see text for details). 



Thus, wc observe a near field regime in which the 
decay of the g 2 appears to be generically more rapid 
than exp(— r/X) followed by the quasi-continuum regime 
where the exponential decay with decay length A is ob- 
served. This demonstrates that the nonaffinity length A, 
indeed, controls the approach toward the expected con- 
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tinuum behavior (i.e., vanishing of the non continuum 1.2 , , , , , — 

L/l = 1 3 9 

modes). Similar results (not presented here for reasons j 29*1 

of space) are also found for the other non-continuum 1 

modes. Specifically, the decay lengths are also found to 

be A within 10%. 8 

Of course, the fixed zero-displacement boundary con- 6 

dition at r — R requires all components of the displace- s° \ 

ment field to vanish there. The essential distinction be- 0.4 - 1 I 

tween the continuum and non-continuum modes lies in * 

■ ■ 

the manner in which their amplitudes decay upon ap- 02 
proach to the rigid boundary. The non-continuum modes 
examined above decay exponentially with a decay length 
proportional to A. We demonstrate in the next section 
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that the continuum mode amplitudes, in contrast, do not 0.2 0.4 0.6 0.8 1 

appear to decay in the same way. In particular, their ap- r/R 
proach to zero at the rigid boundary is not controlled by M 

^ FIG. 7: g ( > versus r/R for X/L » 0.191, R — L and differing 

L/L as shown in the legend. 



2. Continuum modes and finite-size effects 

In this section we address two related points: (i) the 
fundamental change in the structure of the displacement 
field as one moves away from the immediate vicinity of 
the applied force, and (ii) the role of finite size effects in 
our numerical simulation in two-dimensions. We cannot 
address the former point without confronting the latter 
one for the following reason. Due to the presence of the 
rigid boundary, all modes of the displacement field decay 
to zero at the boundary. 

We show in this section, however, that the decay of 
the non-continuum modes is controlled by internal, meso- 
scopic length scales in the network, while the decay of 
the continuum modes is determined by the macroscopic 
geometry of the system including the presence of the 
rigid boundary. Since the distinction between the con- 
tinuum and non-continuum modes depends essentially on 
the presence of the boundary and specifically on the sep- 
aration of the system size R from the internal length 
scales controlling the decay of the non-continuum modes 
(oc A), we consider this distinction along with a more 
general discussion of finite size effects in our simulation. 

To have a well-defined elastic response in two- 
dimensions (meaning that the displacement field vanish 
at large distances from the applied point force) one needs 
to impose a rigid boundary. In order to eliminate the dif- 
ferential influence of the boundary on the various angular 
modes of the displacement field we chose this boundary 
to be a circle of radius R centered on the point of force 
application. This rigid boundary forces all components of 
the displacement field to vanish exactly at r = R. Such a 
rigid boundary can be expected to introduce a boundary 
layer near the edges of our system, which we observe and 
discuss below. 

Figure shows the decay of the <7q mode amplitude 
for systems of size R = L and A « 0.191L. We note that 
the decay of this continuum mode is qualitatively dis- 
tinct from the decay of the non-continuum modes shown 



in Figs.^El andEl We do not observe an exponential de- 
cay of this continuum mode amplitude. Nevertheless, ob- 
serving an exponential decay having a long decay length, 
or, more reasonably, the product of an algebraic function 
and such a weak exponential decay would be difficult to 
resolve in this plot. 

To study these issues further, we plot in Fig.^gp*^ for 
various system sizes from R = L/2 to 5L. We observe 
both an apparent convergence for the larger systems to 
a common curve when distances are measured relative to 
the size of the system, as well as a systematic down-turn 
near the boundary. The first of these observations sug- 
gests that the sample geometry controls the dependence 
of this continuum mode, as opposed to intrinsic lengths 
like A. The down-turn in these data within a distance of 
order A (here, ~ 0.2L) of the boundary is to be expected 
from the observations above concerning the role of the 
length scale A. One can view this crossover length as the 
distance along a filament over which a force applied to 
the filament dissipates/expands into a (quasi-)continuum 
stress/displacement field. Thus, the effects of the bound- 
aries on the individual filaments contacting the boundary 
are expected to propagate at least this distance into the 
system. Because of this, we have also removed all data 
within a distance A of the boundaries in Fig. above, 
which exhibit a similar down-turn near the boundary. 

In Fig. we examine the radial dependence of g£ for 
three different values of X/L for a large system where 
R = 3L. In this case, we see a coincidence of the data 
for different A, indicating that the long-range behavior of 

(r) 

<7o ; is not controlled by A. 

Although we were unable to obtain sufficient statis- 
tics on the non-continuum (and sub-dominant) modes 
for larger systems, the combination of the qualitatively 
different behavior from the continuum modes, together 
with the consistent dependence on the crossover length A 
lead us to conclude that the principal results of the prior 
section are not strongly influenced by finite-size effects. 
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Thus, A controls the disappearance of the non-continuum 
modes of the system. 
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FIG. 8: g ( r) versus r/R for X/L w 0.191, L/l c « 29.1 and 
the system radii i? given in the key. The continuum response 
for the bulk elastic moduli corresponding to these parameter 
values (as given by ©) is shown as a solid line. 
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FIG. 9: (?q versus r/R for R — 3L on log-linear axes, with 
X/L and L/l c as given in the key. Each curve has been shifted 
vertically by a A-dependent scale factor to attempt data col- 
lapse. Since A varies by over a factor of 4, the rough collapse 
is enough to rule out any significant A-dependence on the 
decay. Error bars are not shown for clarity. 

We also see a large discrepancy between the expected 
continuum solution and the observed displacement of the 
simulated hiament network as shown in Fig. 00 This 
difference cannot be simply attributed to finite-size ef- 
fects; as can be seen in Fig. |SJ there is no observable 
convergence of the larger system data for the continuum 
modes to the predictions of continuum elasticity for a 
material having the appropriate Lame coefficients. Nev- 
ertheless, we clearly observe the more rapid decay of the 
non-continuum modes on a length-scale controlled by A. 
One is left with the following puzzle: The tensorial struc- 
ture and the rotational symmetry of each component of 



the displacement field approaches the form required by 
the continuum theory, but the continuum model appears 
to never quantitatively agree with the numerical data. 

We may speculate about the underlying cause of this 
discrepancy. It is clear that that the deformation of the 
material in the immediate vicinity of the point force (this 
"near zone" extends out to a distance oc A from the point 
force) and in a boundary region at the rigid wall is not 
well described by any continuum theory. We suggest that 
under the application of the point force at the origin, the 
network effectively partitions itself into three different 
elastic materials. In the near field region r < A around 
the point force the deformation response to the point 
force is quite complex. Similar complexities appear to 
exist near the rigid wall, where R — r < A . These re- 
gions apply a complex set of tractions on the circles r ~ A 
and r ~ R — r bounding the intermediate region that de- 
forms in a manner consistent with some continuum elas- 
tic material. Because these tractions are not themselves 
determined by a simple, continuum model, the resultant 
deformation of the intermediate region is also not simply 
derivable from an analysis of the elastic Green's function 
for a continuum. 

One might imagine that for very large systems having 
a consequently larger intermediate region, the complex- 
ities of the tractions in the transition zones become less 
significant. Because of the rigid boundary at r — R and 
the low-dimensionality of the system, this may not be 
the case. The continuum Green's function contains log- 
arithmic terms and, due to the rigid boundary, growing 
polynomial terms as well. Thus, we do not expect conver- 
gence to the continuum Green's function even for signifi- 
cantly larger systems. One may ask whether three dimen- 
sional systems will show similarly poor convergence to the 
continuum solutions. Further research here is needed. 

We suggest that we do indeed observe the approach 
of the structure of the elastic Green's function to that 
predicted by continuum theory. We term the region sur- 
rounding the point force where the deformation field has 
the expected form the quasi-continuum. Based on our 
numerical data, we do not expect to find a region in 
which the deformation field agrees quantitatively with 
the predictions of continuum elasticity using the Lame 
coefficients appropriate to the medium as determined by 
uniform stress measurements. This suggests that semi- 
flexible networks admit a highly complex point force re- 
sponse that cannot be fully captured by continuum elas- 
ticity even in the far-field. The full implications of this 
complexity have not been explored. 



3. Energy Density 

Another instructive measure of the monopole response 
is the density of elastic energy p& at a given distance from 
the point perturbation, which measures gradients of the 
displacement field and therefore complements the analy- 
sis of u(x) given above. In addition, since the simulation 
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contains explicit terms for the transverse and longitu- 
dinal filament deformation modes, it is straightforward 
to measure the partitioning of the energy between these 
modes. Based on previous work |5|,|a,|8|, under homo- 
geneous shear strain the partitioning of elastic energy 
between the bending and stretching modes of the fila- 
ments is determined entirely by the affine-to-nonaffine 
crossover. The ratio of L/X was found to control this 
energy partitioning at a macroscopic or average level. It 
remains to be seen how this partitioning of the elastic 
energy occurs in vicinity of a point force. 

The freedom to choose the angle between the point 
force and the direction of the filament to which that force 
is applied allows one to determine locally the partition- 
ing of the elastic energy between bending and stretching 
modes of the filaments. Forces directed along the fila- 
ment axis generate primarily stretching deformations in 
the immediate vicinity of the origin (where the force is 
applied). Forces directed perpendicular to the filament 
axis, however, locally create a large bending deforma- 
tion. As seen in the previous homogenously imposed 
strain deformation calculations, for any given value of 
L/X the network responds very differently to the bend- 
ing or stretching deformations. Thus it is not surprising 
that the decay of the energy density from the point of 
force application to the boundary is strongly dependent 
on whether the force is applied parallel or perpendicular 
to one of the filaments at the crosslink. 

Fig.llOlshows pe for parallel forces, where A is fixed but 
the filament density L/l c varies. There is an approximate 
data collapse onto a single curve as shown. The most 
notable discrepancy occurs near the boundary, where the 
low-density data remains higher than those for larger 
L/l c values. This is most likely due to the boundary layer 
already discussed above for the displacement modes. 
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FIG. 10: Density of elastic energy pE (in arbitrary units) 
versus distance r from an externally applied force that is di- 
rected along the filament to which it is applied. L/l c is varied 
as given in the key, and R = L and X/L ~ 0.191 in all cases. 

In contrast, pe for forces perpendicular to a filament 
exhibits much richer behavior, as shown in Fig. 1111 which 
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FIG. 11: Scaled elastic energy density versus distance from 
the point force for a variety of networks - see legend. In 
each case the force is directed perpendicularly to the filament 
where it is applied. In the upper panel distances are scaled by 
l c and the data has been collapsed in the small r/l c regime. 
The solid line is proportional to exp(— r/l c ), showing / c con- 
trols the initial decay of elastic energy. In the lower panel the 
same energy densities are plotted against distance scaled by 
A. This shows the energy density decays with that character- 
istic length scale at longer ranges. Error bars are not shown 
for reasons of clarity. 



shows the same data plotted against r/l c and r/X. For 
large distances, the data for different L/l c appear to dif- 
fer by only a scale factor. Since there is one arbitrary 
multiplicative factor for each curve (as only the magni- 
tude of the applied force is fixed, not its displacement 
and hence work done on the network), these portions 
of the curves can be made to collapse after scaling, as 
in the parallel force case, suggesting that beyond some 
near-field regime, the decay of elastic energy density is 
once again exponential and governed by A - see Fig. ITT1 
lower panel. However, in this near-field regime for per- 
pendicular forces there is clearly no possibility of such a 
collapse. Instead, it appears that the more rapid decay 
of elastic energy density in this near-field regime is gov- 
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erned by the length l c , as can be seen in the upper panel 
of Fig. ^2 where the (rescaled) energy density data are 
plotted vs. r/l c and shown to decay as ~ exp(— r/l c ) for 
distances r ~ 5l c (ignoring finite size effects). 

Based on the two different regimes of data collapse 
shown in Fig. ^] we conclude that there is a near-field 
region in which the decay of elastic energy is governed 
by a microscopic length scale-the mean distance between 
cross-links-and a longer-range regime in which the spa- 
tial decay of elastic energy density is controlled by the 
mesoscopic length A. To better study this cross-over, we 
examine the partitioning of elastic energy between bend- 
ing and stretching modes of the filaments. Recall that in 
the affine limit, the elastic energy is stored primarily in 
the stretching and compression of the filaments. In the 
nonaffine regime, on the other hand, the elastic energy is 
stored almost entirely in bending modes of the filaments. 

Plotted in Fig. ^] is the proportion of elastic energy 
due to longitudinal filament deformation for the perpen- 
dicular force case. Three data sets are displayed having 
different values of l c and A - see figure caption. The dis- 
tance from the point of force application has been scaled 
by the geometric mean of l c and L leading to the ob- 
served coincidence of the crossover between regimes in 
the data. Close to the point of the application of the 
force, the displacement response is clearly dominated by 
bending modes of the filaments. This is to be expected 
since the perpendicularly directed force directly injects 
bending energy into the system at the origin - the point 
of force application. The energy partitioning is incon- 
sistent with that which is expected based on previous 
homogeneous shear measurements. From that work, one 
expects the network to determine the energy partition- 
ing based only on the ratio: L/X. We observe the frac- 
tion of stretching energy to rapidlyincrease towards the 
expected value based on L/X [51, |8( and we may char- 
acterize each curve has having a 'knee' separating the 
region of varying energy ratio from a nearly constant in- 
termediate regime where the partitioning of elastic ener- 
gies corresponds well with the previously identified frac- 
tion of stretch/compression energy under macroscopic 
strain. This correspondance is demonstrated in Fig. ITS1 
The horizontal lines show the fraction of stretching en- 
ergy in periodic systems subjected to macroscopic shear, 
and clearly coincide with the plateau reached beyond the 
knee. At the largest distances one notices the vanishing 
of bending energy as the fraction of stretching energy 
approaches unity. This last effect is due to the fact that 
each filament has a freely rotating bond at the outer, rigid 
wall. As this wall cannot support torques, the bending 
energy vanishes in a boundary layer whose width is de- 
termined by the mean distance between cross-links. We 
return to this point below. 

The observed coincidence of the knees of all three 
curves under the rescaling of distances by s/l c L suggests 
that this length sets the scale over which injected bending 
energy is redistributed into the combination of bending 
and stretching appropriate for long-length scale deforma- 
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FIG. 12: Proportion of elastic energy due to filament stretch- 
ing, 7^"/(W" + T~L ± ), for three systems having differing values 
of l c (see key) but the same value of A. 
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FIG. 13: The fraction of elastic energy due to filament stretch- 
ing in networks with the same L/l c « 29.1 (in fact, the exact 
same geometries) , R = L but different A. The different values 
of the stretching energy fraction in the intermediate regime 
reflect the differing values of the ratio L/\ for the networks. 
For comparison, the macroscopic energy fraction under bulk 
shear is shown as horizontal dashed lines, in the same order 
as the data points (i.e. X/L pa 0.089, 0.191 and 0.412 from 
top to bottom). 



tions. 

Based on these observations we note that in all net- 
works the strain field acquires the structure of the point 
force response based on continuum elasticity over a (typ- 
ically mesoscopic) length scale of A. When the networks 
is subjected to large, local bending deformations, how- 
ever, it readjusts the partitioning of bending to stretch- 
ing energy over a generally much shorter length scale, l c . 
Thus the system is able to repartition the local clastic 
energy storage to the value appropriate for its L/X ratio 
over smaller length scales than does the system recover 
the expected long length scale structure of its continuum 
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elastic response. 

We have already shown that in the intermediate field 
regime the fraction of stretching elastic energy in the sys- 
tem approaches its far-field value as determined by the 
ratio: L/X. Nevertheless, we in fact find that at the edges 
of our sample the network energy becomes stored solely 
in stretching modes regardless of the value of L/X. As 
mentioned above we attribute this final redistribution of 
the elastic energy density between bending and stretch- 
ing modes to a boundary effect imposed by the freely- 
rotating nature of coupling of the network filaments to 
the rigid boundary. To further test that this final redis- 
tribution is indeed a boundary effect, we consider a few 
larger system sizes as shown in Fig. El In that figure 
the force remains perpendicular to the direction of the 
filament to which it is applied while the system size is 
varied from R = L to R — 5L for a network of fixed 
A. When the data are plotted against the radial distance 
from the point force scaled by system size, we find an 
excellent collapse in the intermediate regime and in the 
putative boundary layer. For the smallest system size 
considered, R = L, we note poorer data collapse in the 
near field region suggesting that one must study systems 
that are are least larger than a single filament length to 
access the bulk behavior of the network with quantita- 
tive accuracy. Clearly, all three curves taken together are 
consistent with the notion of an elastic boundary layer 
that is produced by the freely-rotating boundary at the 
wall and that extends distance approximately equal to A 
into the sample. 



\ i S i I i i i i ■ i I 5 : 



0.1 



0.01 



0.001 



R=L 
R=3L 
R=5L 



0.0001 



0.2 



0.4 0.6 

r/R 



0.( 



FIG. 14: The fraction of elastic energy stored in the H" term 
of the Hamiltonian for perpendicular forces, as a function of 
distance from the origin r for system sizes R = L, 3L and 5L. 
In all cases L/l c w 29.1 and A « 0.19LL. 



B. Network response to force dipoles 

We now consider the mechanical response of the net- 
work to localized force dipoles at the origin. Understand- 
ing this response function is central to elucidating the 



effect of molecular motor activity in the cytoskeleton. 
Many of the general features observed in the response of 
the system to point force monopoles are also in evidence. 
For example in Fig. ll5l we see the rough collapse in the far 
field of the mode amplitude Hq (r) to the continuum so- 
lution for three different networks having the same value 
of A. Note that the amplitude has been scaled by Q since 
the magnitude of the imposed force dipole for each re- 
alization of the network will depend on the distance to 
the next constraint, i.e. cross-link, which is l c . Since the 
displacement field when averaged over network realiza- 
tions is quadrupolar, being the difference in two dipole, 
that length must enter squared. The analogous displace- 
ment field amplitude, h^(r) from the continuum solution 
(solid line in Fig. 115(1 was calculated using unit forces so 
this amplitude is known only up to an overall scale fac- 
tor; that scale factor was adjusted to best fit the data. 
A similar comparison can be made for the other contin- 
uum modes of the displacement field (see the Appendix) . 
Fig. 1161 for example, shows h^lr) with the same arbi- 
trary prefactor. The agreement to the continuum theory 
scaled as discussed above is quite poor. As is also seen in 
the point-force response of the strain field, higher angu- 
lar modes, which vanish in the continuum theory such as 
/i{ are significantly non-zero in the data. The amplitude 
of the non-continuum modes in the dipolar response is 
even more dramatic than in the point-force response of 
the network examined earlier. 

The simulation data for the response of the network 
to force dipoles can be broadly summarized as follows: 
The amplitude of the continuum modes of the dipolar 
displacement field from the naive continuum theory do 
not agree with simulation data. In light of the discussion 
regarding such disagreements between the monopole data 
and the continuum calculation, this is perhaps not sur- 
prising. More interestingly, the observed displacement 
field has significant amplitudes of higher order angular 
modes. In short the universality of the response of semi- 
flexible networks to localized point forces does not appear 
to extend to their response to localized force dipoles. 

We speculate that the principal difference between 
these two cases stems from the fact that the network's 
response to the force dipole probes the more detailed mi- 
croscopic structure of the network in the immediate vicin- 
ity of the point of the force dipole application. The am- 
plitude of the lowest order force multipole communicated 
from the near field region to the far field where we expect 
a continuum based theory to apply is not constrained by 
elementary force balance. Neither are the amplitudes of 
any higher order force multipoles generated within the 
near field region. Thus, upon reaching the inner edge of 
the far field region the force dipole imposed at the origin 
has generated a highly complex set of tractions on the 
rest of the material whose structure depends on local de- 
tails of the connectivity of the network near the origin. In 
contrast for the case of the force monopole, the dominant 
term in those tractions is the fixed total monopolar force 
acting on the intermediate region. The higher order force 
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multipoles created in the inner region decay rapidly with 
distance from the origin leaving one with highly repro- 
ducible results for the response of the network to applied 
forces. The amplitude of the force dipole communicated 
from the near field region to the intermediate region, on 
the other hand, is not similarly constrained. It appears 
that one generically generates large amplitude higher or- 
der force multipoles in addition to whatever force dipole 
is communicated to the intermediate region making con- 
vergence to a simple dipolar form slow and difficult to 
observe in our finite samples. 
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FIG. 15: Hq for dipole forcing for fixed X/L ~ 0.191 compared 
to the continuum solution (smooth line), which has one free 
fitting parameter, namely the overall magnitude of the dipole 
forcing. Both data and curve are negative, so the magnitudes 
have been plotted to allow use of a logarithmic axis. Each 
data curve has been scaled by (L/l c ) 2 to ensure the same 
mean dipole magnitude (see text). The system radius was 
R = L. 



Bulk moduli 



Lastly, we present new results on homogeneous defor- 
mations of semiflexible networks. It has already been 
shown that the macroscopic elastic moduli of this class 
of model networks depend in a crucial way on the ratio 
of A to the hiament length L dH]. For X/L < 1, the 
deformation is approximately affine, whereas non-affine 
deformation modes dominate when X/L ~ 1. Previously 
this was demonstrated only for the shear modulus; here 
we can now confirm that the Young's modulus Y behaves 
in an identical manner. Fig. 1171 shows Y measured from 
uniaxial extension of a rectangular cell, scaled by the pre- 
diction for an affine strain, plotted against X/L for the 
range of L/l c considered in this paper. There is a clear 
data collapse, as for the shear modulus. Furthermore the 
deviation from the affine prediction is small for X/L <C 1, 
but becomes increasingly pronounced as X/L increases. 
This confirms that X/L controls the macroscopic elastic 
response of these systems. 

Fig. 1181 gives the Poisson ratio v for the same systems 
as in Fig. 1171 It is striking to observe that, within error 
bars, v is consistent with the value v = \, which is the 



value expected for an affine deformation |g. However, 
it is apparent from this figure that the measured values 
are consistent with v = \ for all data points, even those 
well into the non-affine regime. The mechanism behind 
this striking robustness currently evades us (is has noth- 
ing to do with incompressibility, which fixes v = 1 in 
two dimensions). Note that v at the rigidity percola- 
tion L/L ^ 5.933 (at which the elastic moduli vanish) 
is ~ g 7], which is clearly inconsistent with the data 
in Fig. 1181 and confirms our earlier claims that the non- 
affine regime is distinct from the scaling regime of the 
transition. 




FIG. 16: h for same systems as in Fig. 1151 The only fitting 
parameter for the continuum solution is the same as used 
previously for hg, so there are no remaining free fitting pa- 
rameters in this plot. Much of the data for r/R > 0.6 is 
negative and hence not visible on these axes. 
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FIG. 17: Young's moduli scaled by the affine prediction versus 
X/L on log-log axes. The affine prediction, which depends 
only on L/l c , can be found in |8(. The symbols are larger 
than the error bars. 
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FIG. 18: The Poisson ratio v for the same systems as in 
Fig. ^3 The solid line corresponds to v = ~. 



IV. DISCUSSION 

In this paper we have presented the results of numerical 
studies on the response of semiflexible networks to both 
point forces and to homogeneously imposed strain. The 
data presented on the Young's modulus taken in combi- 
nation with previous work on the static shear modulus 
shows that the mechanical response of the system can be 
understood in terms of Lame coefficients that depend on 
the ratio of the filament length to the nonaffinity length: 
L/X. The Poisson ratio of the material, however, appears 
to be remarkably insensitive to this ratio. We can offer 
no explanation for this insensitivity at this time. The 
data presented on the point force response form a neces- 
sary compliment to previous work on the development of 
a long length scale elastic theory of such materials. 

Based on these data, it appears that the storage of elas- 
tic energy and the structure of the strain field is rather 
complex in the immediate vicinity of the applied point 
force. We characterize these quantities by considering 
three qualitatively different regimes as a function of ra- 
dial distance from the applied point force. Immediately 
surrounding the point force in the near-field regime we 
find that the partitioning of elastic energy into bending 
and stretching is determined primarily by the angle be- 
tween the applied force and the direction of the filament 
to which that force is applied. The disorder-averaged 
strain field is rather complex having higher angular har- 
monics that predicted by continuum elasticity theory. In 
the intermediate field regime farther from the point force 
the partitioning of elastic energy is determined solely by 
the ratio L/X as found in the homogeneous shear mea- 
surements. The higher angular harmonics present in the 
strain field appear to decay exponentially with a decay 
length proportional to A with a constant of proportional- 
ity near unity. These results taken in combination with 
our previous work suggests that one may think of A as 
setting the minimum length scale for the applicability of 



continuum modeling quite generally. 

In the quasi-continuum regime (r > A), we find a strain 
field consistent with structure of that predicted by con- 
tinuum model. The strain itself, however, cannot be sim- 
ply computed from a knowledge of the effective Lame co- 
efficients. We believe the source of this discrepancy is the 
fact that the intermediate region, being poorly described 
by the continuum theory, applies a much more complex 
set of tractions to the material in the far-field where the 
continuum theory must apply. If these boundary con- 
ditions were known, we suspect that one could in fact 
calculate the resulting displacement field using the con- 
tinuum theory. This belief is supported by the fact that 
in the far-field regime the displacement field for different 
networks having the same A collapse onto the same curve 
as one would expect for any effective continuum theory 
in which the Lame coefficients depend on A. 

To summarize these results, we suggest that there is 
an emerging description of the mechanics semiflexible 
networks. There mechanical behavior over length scales 
longer than A appears to be described by a modified elas- 
tic theory in which the effective Lame constants depend 
on the ratio L/X. For point forces and presumably any 
force applied over regions having a characteristic length 
scale that is small compared to A, the local response of 
the network is quite complex and the material appears to 
be anomalously compliant; at longer length scales, how- 
ever, the structure of the deformation field appears to 
be consistent with the predictions of continuum elastic- 
ity. We believe that it should be possible to construct an 
elastic continuum theory of these networks that is appli- 
cable in both the intermediate and far-field and that is 
based on a modified gradient expansion of the strain field 
incorporating explicitly the mesoscopic length A. The 
behavior of the strain field on scales much smaller than 
A appears to depend on other, more microscopic length 
scales. 

Understanding the response of semiflexible networks 
to localized forces is a necessary for both microrheolog- 
ical investigations of semiflexible networks and for un- 
derstanding the effect of molecular motors in the cy- 
toskeleton. Clearly, further numerical investigations are 
required as well as a theoretical examination of the de- 
velopment of elastic continuum models applicable to the 
intermediate and far-field regimes. 

Moreover, additional investigations are required to ex- 
amine the analogous questions in three dimensional semi- 
flexible networks. While we expect the basic physics 
outlined above, including the existence of a mesoscopic 
length A, to persist in three dimensions, the expected 1/r 
decay of the elastic Green's function should alter the re- 
sults. Moreover the more rapid decay of the displacement 
field and energy density in the continuum three dimen- 
sional system may further simplify the structure of the 
analogous Green's functions for the semiflexible network. 
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APPENDIX 

Here we derive the displacement field u(x) predicted 
by continuum elasticity for the cases of both monopole 
and dipole forcing. For the monopole case, a point force 
f<5(x) is applied to the origin of an elastic sheet with Lame 
coefficients /i and A, or equivalently [i and the Poisson 
ratio v = A/(2/x + A). For an isotropic elastic body, the 
stress obeys — \i(diUj + djui) + XSg V • u and the 
resulting equation for force balance is 



+ A)9i(V-u)+/iV 2 



/i*(x) 



(7) 



Solving this in polar coordinates (r, <fi) with the boundary 
condition u = at a radius R from the origin eventually 
leads to 

< on ° = ^-{ifMr/R)-^ifjfi 
-i 



c 2 [r/RY 
' (r/RY 



Zrirjfj - c 2 fi 



(8) 



where c\ = (1 + v) /(3 — v) and c 2 = (5 + f)/(3 — v). Then 
the only two non-zero modes according to the definition 
of the <3>W given in Q are 



^l-Ar, +2c 2 (r/R)- 2 -2(r/R) 2 }, (9) 

9 { f) = ^{±Hr/R)~c 2 {r/R)- 2 + c 2 (r/R) 2 ). 

(10) 

A naive calculation of the corresponding dipole solu- 
tion would simply superpose the above monopole solution 
for two point forces, f<5(x) and — f<S(x — e) with e = ef , 
and then take the limit e — > 0. However, this ignores the 
boundary at radius _R, which should be kept fixed but is 
shifted a distance e by the above procedure. An exact 
calculation would require the monopole solution for force 
applied near to (but not at) the center of a circular sys- 
tem, which, since it no longer obeys radial symmetry, is 



likely to be highly complex. Here we ignore such issues 
and simply use the above monopole solution, in the ex- 
pectation that it will closely approximate the exact case 
except possibly near the boundary. The displacement 
field induced by the force dipole is then 



dip 

< 1 — 



[4ci - 2c 2 {r/R)- 2 + 2(r/R) 2 } 

+ 8fi(? -f) 2 [c 2 (r/R)- 2 - Cl ] 

+ 2/i(f -f) [2( Cl -l)-2c 2 (r/i?)- 2 

-(c 2 -l)(r/R) 2 ]\. (11) 



As explained in Sec. Ill CI the dipole solution above ap- 
plies for a single realization on large length scales, but 
after averaging over many networks on short lengths (as 
in the simulations) the dipole moment vanishes, leaving a 
quadrupole displacement field. The required quadrupole 
is one consisting of two parallel dipoles of equal and op- 
posite magnitude, aligned along their axes, i.e. 



quad 



-V 

e 



dip 



Stt/i r 2 



j/ 4 [4(2ci - 1) - 6c 2 (r/R)- 2 

+ 2(2-c 2 )(rAR) 2 ] 
+ 24f 4 (f -f) [c 2 (r/R)- 2 - Cl ] 
+ 8/j(r -f) 2 [l-2 Cl +3c 2 (r/i?)- 2 ] 

+ 16fi(r • f) 3 [2 Cl - 3c 2 (r/i?)- 2 ] | . (12) 

This gives the displacement field in response to known 
forces of magnitude /. Since it is rather the displacement 
that is controlled, / is a free parameter. Finally, the non- 
zero modes are 



,(/) 



,(/) 



Maffinc 


1 


TTfl 




Maffinc 


1 


27T/i 




Maffinc 


1 


47T/Z 




Maffine 


1 


47T/X 





{-ci}, (13) 

{2 Cl -3c 2 (r/i?)- 2 }, (14) 
{3c 2 (r/i?)- 2 

+ (2-c 2 )(r/ J R) 2 }, (15) 

3c 2 (r/R)- 2 }. (16) 
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